Notebook 3: Join Tables and Make Maps

Introduction

In this exercise we’ll practice joining climate data tables from Cal-Adapt to other tables so we can create indices and map results. Our exercise will focus on census tracts in Kern County, and the end goal will be to create maps with a composite index of census tract cumulative impact score from CalEnviroScreen (a measure of exposure to pollution) and the anticipatd increase in temperature from historic condition to the end-of-century under RCP8.5.

Setup

First we load the packages we need for this exercise, which include caladaptr, dplyr, and sf:

library(caladaptr)
library(sf)
library(dplyr)
library(lubridate)
library(units)
library(ggplot2)
library(tidyr)
library(scales)
library(tmap)

Get All the Census Tracts for Kern County

Start by downloading all the census tracts:

tracts_sf <- ca_aoipreset_geom("censustracts")

Reading layer `censustracts' from data source `C:\Users\Andy\AppData\Local\R\cache\R\caladaptr\censustracts.gpkg' using driver `GPKG'
Simple feature collection with 8034 features and 63 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -13850030 ymin: 3833633 xmax: -12705030 ymax: 5162403
projected CRS:  WGS 84 / Pseudo-Mercator

head(tracts_sf)

Simple feature collection with 6 features and 63 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -13610980 ymin: 4556079 xmax: -13604620 ymax: 4563255
projected CRS:  WGS 84 / Pseudo-Mercator
       tract pop2010 ciscore ciscorep ces_3_0_percentile_range ozone ozonep pm pmp
1 6001400100    2937       3        2     1-5% (lowest scores)     0      8  8  31
2 6001400200    1974       1        0     1-5% (lowest scores)     0      8  8  31
3 6001400300    4865      12       18                   15-20%     0      8  8  31
4 6001400400    3703      10       15                   15-20%     0      8  8  31
5 6001400500    3517      18       35                   30-35%     0      8  8  31
6 6001400600    1571      18       34                   30-35%     0      8  8  31
  diesel dieselp drink drinkp pest pestp rseihaz rseihazp traffic trafficp cleanups
1     27      81    70      4    0     0     485       50     929       64        4
2     42      94    70      4    0     0     442       49    1392       81        0
3     42      94    70      4    0     0     426       48    1207       75        0
4     42      94    70      4    0     0     444       49    1153       74        0
5     42      94    70      4    0     0     448       49     668       48        3
6     42      94    70      4    0     0     426       48    1743       88        1
  cleanupsp gwthreats gwthreatsp haz hazp iwb iwbp swis swisp pollution pollutionscore
1        42         4         27   6   98   2   29    1    33        37              4
2         0        20         72   0    0   0    0    0     0        30              3
3        13        35         85   0    0   0    0    0     0        31              3
4         0        38         87   0    0   0    0    0     0        30              3
5        35        21         73   0    0   0    0    0     0        29              3
6        24        20         72   0    0   0    0    0     0        32              3
  pollutionp asthma asthmap lbw lbwp cvd cvdp edu edup ling lingp pov povp unemp unempp
1         40     18       6   2    9   3    1   1    4    3    29   7    4    NA     NA
2         20     24      15   1    0   3    5   1    3    0     0  11    9   2.2      1
3         23     37      38   4   29   4    8   5   20    8    55  14   15   8.8     46
4         21     60      70   3   26   4   12   6   24    2    17  16   20   3.3      4
5         18    110      95   5   68   7   44   3   11    2    17  25   36   6.5     25
6         25    122      97   1    1   8   55   7   31    1     8  25   37  14.9     84
  housingb housingbp popchar popcharscore popcharp children_under10_pct
1        8         8       8            0        1                    9
2        4         1       4            0        0                   12
3       16        40      30            3       20                   11
4       15        36      27            2       17                   11
5       19        59      49            5       49                    9
6       13        26      44            4       41                    7
  pop_11_64_years_pct elderly_over_65_pct hispanic_pct white_pct african_american_pct
1                  69                  21            4        70                    4
2                  71                  16            7        78                    1
3                  78                  10            8        66                   10
4                  79                   9            9        65                   12
5                  82                   8            9        50                   26
6                  82                   9            8        42                   39
  native_american_pct asian_american_pct other_pct cidecile civigintile id
1                   0                 15         4        1           1  1
2                   0                  7         5        1           1  2
3                   0                  8         5        2           4  3
4                   0                  7         5        2           4  4
5                   0                  6         6        4           7  5
6                   0                  5         4        4           7  6
                            geom
1 MULTIPOLYGON (((-13608465 4...
2 MULTIPOLYGON (((-13609690 4...
3 MULTIPOLYGON (((-13610548 4...
4 MULTIPOLYGON (((-13610124 4...
5 MULTIPOLYGON (((-13610979 4...
6 MULTIPOLYGON (((-13610819 4...

We can pull out just the tracts for Kern County by decomposing the tract id. All census tracts in Kern County will have an idea number that starts with “6029” (06 is for California, 029 is Kern County’s FIPs code).

kern_tracts_sf <- tracts_sf %>%
  filter(tract %>% as.character() %>%  stringr::str_starts("6029"))

View the attribute table. You’ll notice the census tracts have values from CalEnviroScreen, including the cumulative impact score (ciscore), which we’ll come back to you later.

head(kern_tracts_sf)

Simple feature collection with 6 features and 63 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -13260680 ymin: 4217755 xmax: -13248100 ymax: 4228472
projected CRS:  WGS 84 / Pseudo-Mercator
       tract pop2010 ciscore ciscorep ces_3_0_percentile_range ozone ozonep pm pmp
1 6029000101   12240      47       86                   85-90%     0     91 18 100
2 6029000102    3055      33       66                   65-70%     0     91 18 100
3 6029000200    7644      50       89                   85-90%     0     91 18 100
4 6029000300    4217      47       86                   85-90%     0     91 18 100
5 6029000400    4319      57       95 95-100% (highest scores)     0     91 18 100
6 6029000503    6347      15       28                   25-30%     0     91 18 100
  diesel dieselp drink drinkp pest pestp rseihaz rseihazp traffic trafficp cleanups
1     21      66   485     52  137    78      87       24     512       32       31
2     21      65   288     31    0     0      21       12     209        6       11
3     21      67   288     31    0     0      53       20     528       34       21
4     21      66   288     31    0     0      49       19     406       22       17
5     49      96   786     89    0     0      89       24     761       54       14
6     26      79   691     82 1352    91     496       51     984       67        1
  cleanupsp gwthreats gwthreatsp haz hazp iwb iwbp swis swisp pollution pollutionscore
1        91        80         96  11  100   0    0    2    50        64              7
2        65        29         81   0   43   0    0    0     0        41              5
3        85        24         77   3   95   0    0    0     0        49              6
4        81         8         43   2   93   0    0    0     0        45              5
5        75         4         29   2   93   0    0    0     0        56              6
6        18         0          4   0   26   0    0    0    20        57              7
  pollutionp asthma asthmap lbw lbwp cvd cvdp edu edup ling lingp pov povp unemp unempp
1         96     80      86   5   61  11   87  11   44    1    10  30   46   8.7     45
2         50     80      86   5   64  11   87  15   53    0     0  49   72  17.4     91
3         72     80      86   6   86  11   87  22   67    2    14  75   97  16.7     90
4         60     80      86   7   96  11   87  26   71    3    23  68   92  23.6     98
5         86     79      85   6   80  11   84  24   70    4    35  76   97  21.3     97
6         88     33      31   2    4   7   46  10   42    1    12  11    9   5.3     15
  housingb housingbp popchar popcharscore popcharp children_under10_pct
1       16        44      57            6       62                   17
2       12        24      63            6       70                   15
3       33        93      79            8       92                   16
4       24        74      80            8       93                   16
5       32        91      80            8       93                   17
6        6         4      21            2       10                   14
  pop_11_64_years_pct elderly_over_65_pct hispanic_pct white_pct african_american_pct
1                  74                   7           19        74                    0
2                  69                  15           11        82                    0
3                  75                   7           21        73                    0
4                  71                  12           21        73                    0
5                  73                   9           21        72                    1
6                  76                   8           17        74                    0
  native_american_pct asian_american_pct other_pct cidecile civigintile  id
1                   1                  2         2        9          18 965
2                   1                  0         3        7          14 966
3                   1                  0         3        9          18 967
4                   1                  0         2        9          18 968
5                   1                  0         3       10          20 969
6                   0                  3         2        3           6 970
                            geom
1 MULTIPOLYGON (((-13256411 4...
2 MULTIPOLYGON (((-13249411 4...
3 MULTIPOLYGON (((-13253779 4...
4 MULTIPOLYGON (((-13252464 4...
5 MULTIPOLYGON (((-13253750 4...
6 MULTIPOLYGON (((-13260684 4...

For mapping, we’ll also import the county boundary for Kern:

kern_bnd_sf <- ca_aoipreset_geom("counties") %>% 
  filter(fips == "06029")

Reading layer `counties' from data source `C:\Users\Andy\AppData\Local\R\cache\R\caladaptr\counties.gpkg' using driver `GPKG'
Simple feature collection with 87 features and 54 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -13871160 ymin: 3833648 xmax: -12625080 ymax: 5416187
projected CRS:  WGS 84 / Pseudo-Mercator

Plot them to make sure everything looks right:

library(tmap)
tmap_mode("plot")
tm_shape(kern_tracts_sf) +
  tm_polygons() +
tm_shape(kern_bnd_sf) +
  tm_borders(col = "red", lty = 2) +
tm_layout(main.title = "Census Tracts in Kern County, CA")

Now we can create two API requests for the historic values and end of century values:

## Create an API request for historic values of tasmax, using a 32 model ensemble
kern_trcts_hist_cap <- ca_loc_aoipreset(type = "censustracts",
                                   idfld = "tract",
                                   idval = kern_tracts_sf %>% pull(tract)) %>%
  ca_cvar("tasmax") %>%
  ca_period("year") %>%
  ca_gcm("ens32max") %>%
  ca_scenario("historical") %>%
  ca_years(start = 1985, end = 2005) %>%
  ca_options(spatial_ag = "mean")

kern_trcts_hist_cap

Cal-Adapt API Request
Location(s): 
  AOI Preset: censustracts
  tract(s): 6029000101, 6029000102, 6029000200, ...
Variable(s): tasmax
Temporal aggregration period(s): year
GCM(s): ens32max
Scenario(s): historical
Slug(s): NA
Dates: 1985-01-01 to 2005-12-31
Options: 
  spatial ag: mean
  temporal ag (add'l): NA
 

plot(kern_trcts_hist_cap)

Fetch the historic data:

kern_trcts_hist_tbl <- kern_trcts_hist_cap %>% ca_getvals_tbl(quiet = TRUE)

## backup: kern_trcts_hist_tbl <- readRDS("data/kern_trcts_hist_tbl.rds")

dim(kern_trcts_hist_tbl)

[1] 3171    8

head(kern_trcts_hist_tbl)

Next compute the average temperature for each tract.

kern_trcts_mean_temp_hist <- kern_trcts_hist_tbl %>%
  group_by(tract) %>%
  summarise(mean_temp_hist = mean(val))

`summarise()` ungrouping output (override with `.groups` argument)

kern_trcts_mean_temp_hist

Do the same for the end-of-century period with RCP85.

kern_trcts_prj_cap <- ca_loc_aoipreset(type = "censustracts",
                                        idfld = "tract",
                                        idval = kern_tracts_sf %>% pull(tract)) %>%
  ca_cvar("tasmax") %>%
  ca_period("year") %>%
  ca_gcm("ens32max") %>%
  ca_scenario("rcp85") %>%
  ca_years(start = 2070, end = 2099) %>%
  ca_options(spatial_ag = "mean")

kern_trcts_prj_cap

Cal-Adapt API Request
Location(s): 
  AOI Preset: censustracts
  tract(s): 6029000101, 6029000102, 6029000200, ...
Variable(s): tasmax
Temporal aggregration period(s): year
GCM(s): ens32max
Scenario(s): rcp85
Slug(s): NA
Dates: 2070-01-01 to 2099-12-31
Options: 
  spatial ag: mean
  temporal ag (add'l): NA
 

Fetch values:

kern_trcts_prj_tbl <- kern_trcts_prj_cap %>% ca_getvals_tbl(quiet = TRUE)

## backup: kern_trcts_prj_tbl <- readRDS("data/kern_trcts_prj_tbl.rds")

dim(kern_trcts_prj_tbl)

[1] 4530    8

head(kern_trcts_prj_tbl)

Compute the mean for each tract:

kern_trcts_mean_temp_prj <- kern_trcts_prj_tbl %>%
  group_by(tract) %>%
  summarise(mean_temp_prj = mean(val))

`summarise()` ungrouping output (override with `.groups` argument)

Now we can join the tables:

dim(kern_tracts_sf)

[1] 151  64

dim(kern_trcts_mean_temp_hist)

[1] 151   2

dim(kern_trcts_mean_temp_prj)

[1] 151   2

kern_tracts_plus_temps_sf <-
  kern_tracts_sf %>%
  left_join(kern_trcts_mean_temp_hist, by = "tract") %>%
  left_join(kern_trcts_mean_temp_prj, by = "tract") %>%
  mutate(temp_increase = mean_temp_prj - mean_temp_hist) %>%
  select(tract, ciscore, temp_increase)

head(kern_tracts_plus_temps_sf)

Simple feature collection with 6 features and 3 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -13260680 ymin: 4217755 xmax: -13248100 ymax: 4228472
projected CRS:  WGS 84 / Pseudo-Mercator
       tract ciscore temp_increase                           geom
1 6029000101      47  5.023285 [K] MULTIPOLYGON (((-13256411 4...
2 6029000102      33  5.034754 [K] MULTIPOLYGON (((-13249411 4...
3 6029000200      50  5.016184 [K] MULTIPOLYGON (((-13253779 4...
4 6029000300      47  5.016184 [K] MULTIPOLYGON (((-13252464 4...
5 6029000400      57  5.016184 [K] MULTIPOLYGON (((-13253750 4...
6 6029000503      15  5.011804 [K] MULTIPOLYGON (((-13260684 4...

Let’s make a choropleth map of the Cumulative Impact score as well as the mean temp_increase:

library(tmap)
tmap_mode("plot")

tmap mode set to plotting

tm_shape(kern_tracts_plus_temps_sf) +
  tm_polygons (col = "ciscore",
               n = 10,
               style = "cont",
               palette = "YlOrRd",
               colorNA = "grey50",
               legend.reverse = TRUE,
               title = "CI Score"
  ) +
  tm_layout(main.title = "Kern County Cummulative Impact Score\nCalEnviroScreen",
            main.title.size = 0.9,
            legend.position = c("left", "bottom")) +
  tm_scale_bar(position = c("right", "bottom"))

Mean temp increase:

tm_shape(kern_tracts_plus_temps_sf) +
  tm_polygons (col = "temp_increase",
               n = 10,
               style = "cont",
               palette = "YlOrRd",
               colorNA = "grey50",
               legend.reverse = TRUE,
               title = "Mean Temp Increase (K)"
  ) +
  tm_layout(main.title = "Kern County Mean Temp Increases\nHistoric Period - End of Century, 32-ens GCM, RCP85",
            main.title.size = 0.9,
            legend.position = c("left", "bottom")) +
  tm_scale_bar(position = c("right", "bottom"))

NA
NA

To plot CI score and mean temp increase together, we can rescale each one 0..1 and then multiply them together.

kern_tracts_plus_temps_idx_sf <-
  kern_tracts_plus_temps_sf %>%
  mutate(ciscore_01 = scales::rescale(ciscore),
         temp_increase_01 = scales::rescale(as.numeric(temp_increase))) %>%
  mutate(csi_temp_idx = ciscore_01 * temp_increase_01)

head(kern_tracts_plus_temps_idx_sf)

Simple feature collection with 6 features and 6 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -13260680 ymin: 4217755 xmax: -13248100 ymax: 4228472
projected CRS:  WGS 84 / Pseudo-Mercator
       tract ciscore temp_increase                           geom ciscore_01
1 6029000101      47  5.023285 [K] MULTIPOLYGON (((-13256411 4...  0.6025641
2 6029000102      33  5.034754 [K] MULTIPOLYGON (((-13249411 4...  0.4230769
3 6029000200      50  5.016184 [K] MULTIPOLYGON (((-13253779 4...  0.6410256
4 6029000300      47  5.016184 [K] MULTIPOLYGON (((-13252464 4...  0.6025641
5 6029000400      57  5.016184 [K] MULTIPOLYGON (((-13253750 4...  0.7307692
6 6029000503      15  5.011804 [K] MULTIPOLYGON (((-13260684 4...  0.1923077
  temp_increase_01 csi_temp_idx
1        0.1280552   0.07716145
2        0.1439292   0.06089311
3        0.1182271   0.07578663
4        0.1182271   0.07123943
5        0.1182271   0.08639676
6        0.1121651   0.02157021

Plot our index:

tm_shape(kern_tracts_plus_temps_idx_sf) +
  tm_polygons (col = "csi_temp_idx",
               n = 10,
               style = "cont",
               palette = "YlOrRd",
               colorNA = "grey50",
               legend.reverse = TRUE,
               title = "CI Score x Temp Increase"
  ) +
  tm_layout(main.title = "Kern County Cummulative Impact Score * End-of-Century Temp Increase \nCalEnviroScreen",
            main.title.size = 0.9,
            legend.position = c("left", "bottom")) +
  tm_scale_bar(position = c("right", "bottom"))

Conclusion

In this example we saw how to join data from Cal-Adapt to other tables using the feature id column. This gives us the ability to combine different types of location-based information for analysis and visualization.